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A  detailed  collisional-radiative  model  is  developed  and  coupled  with  a  single-fluid,  two-temperature  convec¬ 
tion  model  for  the  transport  of  shock-heated  argon.  The  model  is  used  in  a  systematic  approach  to  examine 
the  effects  of  the  collision  cross  sections  on  the  shock  structure,  including  the  relaxation  layer  and  subse¬ 
quent  radiative-cooling  regime.  We  present  a  comparison  with  previous  experimental  results  obtained  at  the 
University  of  Toronto  and  the  Australian  National  University,  which  serve  as  benchmarks  to  the  model.  It 
is  shown  that  ionization  proceeds  via  the  ladder-climbing  mechanism  and  is  dominant  from  the  upper  levels 
as  compared  to  the  metastable  states.  Taking  this  into  account,  the  present  model  is  able  to  accurately 
reproduce  the  metastable  populations  in  the  relaxation  zone  measured  in  previous  experiments,  which  is  not 
possible  with  a  two-step  model.  Numerical  results  of  the  radiative-cooling  region  are  in  close  agreement  with 
experiments  and  have  been  obtained  without  having  to  consider  radiative  transport.  In  particular,  sponta¬ 
neous  free-bound  emission  to  the  upper  levels  together  with  Bremsstrahlung  emission  account  for  nearly  all 
radiative  losses;  all  other  significant  radiative  processes,  resulting  in  a  transitions  into  the  ground-state,  are 
mostly  self-absorbed  and  have  a  lesser  impact.  The  effects  of  electron  heat  conduction  are  also  considered  and 
shown  to  have  a  large  impact  on  the  electron-priming  region  immediately  behind  the  shock  front;  however, 
the  overall  effect  on  induction  length  is  small. 

Keywords:  Non-Equilibrium  Collisional-Radiative  Plasma  Shocks 


I.  INTRODUCTION 

Non-equilibrium  effects  in  shock-heated  plasma  pro¬ 
vide  valuable  information  such  as  radiative  heat  fluxes 
and  associated  spectra  and  are  consequently  of  great  in¬ 
terest  in  research  areas  such  as  hypersonic  aerodynamics 
and  re-entry  flows,  as  well  as  other  high-enthalpy  plasma 
dynamic  applications.  Without  an  energy  sink  for  rota¬ 
tional  or  vibrational  modes  or  dissociative  processes,  high 
shock  speeds  and  temperature  jumps  are  easily  obtained 
with  noble  gases,  making  it  ideal  to  study  electronic  ex¬ 
citation  and  ionization  kinetics  towards  the  validation  of 
collisional-radiative  (CR)  models.  The  study  of  multi¬ 
dimensional  ionizing  shock  structures  is  a  first  step  to¬ 
wards  our  ultimate  objective  of  a  fully-coupled,  multi¬ 
dimensional  flow  code  for  the  analysis  of  complex  non¬ 
equilibrium  plasma,  which  allows  us  to  test  high-order 
transport  solvers  to  resolve  the  dynamical  nature  of  the 
shock  structure  as  well  as  the  detailed  kinetics  required 
to  solve  the  non-equilibrium  processes  in  the  relaxation 
region.  While  the  relaxation  of  shock-heated  Argon  has 
been  the  focus  of  many  numerical  investigations,  the  CR 
kinetics  employed  have  been  modest  at  best.  The  first 
experiments  on  ionization  relaxation  in  argon  were  car¬ 
ried  out  by  Petschek  and  Byron1,  the  results  of  which 
were  used  to  benchmark  theoretical  models  in  subsequent 
work,  while  the  structure  of  an  ionizing  shock  wave  in 
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argon  was  first  calculated  by  Biberman  and  Yakubov2. 
The  majority  of  computations  have  been  based  on  the 
two-step  reaction  mechanism  proposed  by  Weymann3  in 
which  Argon  atoms  are  initially  excited  to  the  lowest  ex¬ 
cited  levels  by  atom  impact  from  which  the  valence  elec¬ 
tron  is  ejected  by  further  impact,  while  Wojciechowski 
and  Weymann4  showed  that  the  “ladder-climbing”  ion¬ 
ization  proposed  by  de  Boer  and  Grimwood5  becomes  in¬ 
creasingly  important  as  gas  temperatures  and  pressures 
increase.  Coupling  of  the  excitation  and  ionization  ki¬ 
netics  with  the  fluid  dynamics  is  critical  to  an  accurate 
analysis  of  the  plasma,  and  yields  further  complexity  to 
the  problem.  For  example,  McLaren  and  Hobson6  took 
into  account  the  effect  of  the  bounday  layer  on  the  over¬ 
all  relaxation  length  according  to  Mirel’s  theory7,8.  If  the 
boundary  layer  is  small  relative  to  the  cross  section  of  the 
shock  tube,  it  is  usually  argued  that  the  shock  structure 
can  be  well  approximated  as  a  one-dimensional  structure. 
Nevertheless,  we  will  show  in  a  companion  paper  that 
multi-dimensional  effects  can  be  important.  Similarly, 
unsteady  dynamical  effects  can  result  from  kinetics  and 
fluid  coupling:  shock  instabilities  in  argon  were  noted  as 
far  back  as  Shreffler  &  Christian9,  and  a  previous  numer¬ 
ical  study10  demonstrated  the  effect  of  an  unsteady  wave 
action  within  the  induction  zone,  and  suggested  a  rele¬ 
vance  to  the  experimental  observations  of  Glass  &  Liu11. 
The  current  work  expands  on  that  study  by  using  a  more 
accurate  and  detailed  kinetics  model,  and  demonstrating 
multi-dimensional  and  unsteady  effects  in  a  companion 
article. 

Below,  the  focus  is  placed  on  the  collisional-radiative 
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TABLE  I.  Atom-atom  cross  section  parameters  for  excitation 
from  ground  state  used  in  various  attempts  to  reproduce  the 
shock  structure  in  argon. 


model  itself  and  the  timing  (calibration)  of  some  of  the 
rates  and  cross-sections,  using  features  of  the  experi¬ 
mental  data  from  ionizing  shocks  such  as  the  relaxation 
length  and  ionization  fraction.  In  particular,  it  has  been 
found  that  atom-atom  cross  sections  for  excitation  from 
ground  state  have  the  greatest  influence  on  the  shock 
structure  and  the  relaxation  length.  Since  the  gas  im¬ 
mediately  downstream  of  the  shock  is  not  yet  ionized, 
atom-atom  collisions  behind  the  shock  discontinuity  gen¬ 
erate  the  initial  priming  electrons,  which  in  turn  become 
the  dominant  ionizing  species  once  sufficient  number  den¬ 
sities  have  been  obtained.  Some  values  of  the  atom-atom 
cross  sections  used  in  previous  studies,  listed  in  Table  I, 
show  that  the  range  of  values  exceeds  an  order  of  mag¬ 
nitude  and  indicates  the  relative  uncertainty  in  this  par¬ 
ticular  cross  section. 

We  point  out  that  with  the  exception  of  the  UTIAS 
experiments,  the  shock-tubes  were  relatively  small,  the 
largest  having  an  inner  diameter  of  5cm.  In  such  shock 
tubes,  the  viscous  boundary  layer  becomes  a  significant 
portion  of  the  flow,  consequently  reducing  the  overall  re¬ 
laxation  length.  McLaren  and  Hobson6  took  this  into 
account  using  Mirel’s  transformation  and  found  signifi¬ 
cant  reduction  in  the  atom-impact  excitation  cross  sec¬ 
tions  in  comparison  with  previous  work.  Based  on  the 
UTIAS  data,  Kaniel17  showed  that  the  slope  of  the  elec¬ 
tron  avalanche  was  better  reproduced  by  increasing  the 
electron-impact  cross  sections.  This  increase  in  electron- 
impact  excitation  cross  sections  should  be  accompanied 
by  a  decrease  in  the  atom-impact  cross  sections.  In  this 
first  article  we  compute  steady-state  shock  structures 
and  comparison  with  experiments  shows  some  evidence 
that  the  ionizing  shock  structure  can  be  used  to  place 
bounds  on  the  atom-impact  cross-sections,  while  further 
constraints  can  be  obtained  using  unsteady  effects  (dis¬ 
cussed  in  companion  article). 


II.  GOVERNING  EQUATIONS 

For  the  conditions  of  interest,  the  plasma  is  sufficiently 
collisional  that  a  single-fluid  continuum  model  can  be  as¬ 


sumed.  Thermal  non-equilibrium  effects  are  accounted 
for  by  a  multi-temperature  (2T)  formulation,  while  the 
electronic  states  are  transported  as  individual  species,  al¬ 
lowing  a  non-Boltzmann  description  of  the  Atomic  State 
Distribution  Function  (ASDF).  The  convective  trans¬ 
port  terms  form  a  hyperbolic  system  of  equations  which 
are  written  in  conservative  form  and  solved  with  a  3rd- 
order  finite- volume,  monotonic  numerical  scheme.  In 
this  article,  only  steady-state  solutions  are  presented, 
but  the  same  multi-fluid  equations  apply  to  the  multi¬ 
dimensional  and  unsteady  calculations  presented  in  a  sec¬ 
ond  article.  The  differential  form  of  the  governing  equa¬ 
tions  is,  in  the  absence  of  applied  fields  (here  in  ID): 


dQ 

dt 


(1) 


where  Cl  describes  all  the  collisional-radiative  kinetics  (as 
well  as  energy  exchange  from  elastic  collisions)  while  Q 
and  F  are  the  vectors  of  conserved  variables  and  fluxes, 
respectively. 
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being  the  total  energy  density  of  the  plasma  and  H  =  E+ 
p  the  total  enthalpy  density.  In  Eq.  3,  £h  =  T,s^epses/ p 
is  the  total,  specific  internal  energy  density  of  the  heavy 
fluid.  Note  also  that  Eq.  2  contains  all  3  components 
of  the  velocity  (u,v,w)  although  only  one  component  is 
non-zero  in  the  ID  geometry  studied  here.  The  electron 
energy  equation  has  been  expressed  in  terms  of  the  en¬ 
tropy  function  se  =  pe/ ple  (cf.18)  in  order  to  maintain  a 
strict  divergence  form  for  the  convected  quantities.  In¬ 
troducing  the  flux  Jacobian  A  =  dF/dQ ,  Eq.  1  can  be 
written  in  quasilinear  form: 


dQ 

dt 


,dQ 

dx 


=  Cl 


(4) 


In  the  absence  of  discontinuities  (i.e.  behind  the  shock), 
this  form  will  be  used  for  steady-state  solutions,  as  de¬ 
scribed  below.  Closure  of  the  system  is  obtained  with 
an  equation  of  state  (EOS):  the  plasma  is  assumed  here 
to  be  ideal  and  thermally  perfect,  together  with  Dalton’s 
law  of  partial  pressures,  i.e.  p  =  nfcsT  and  es  =  f(T ), 
such  that  des  =  Cv^s{T)dT  and  CVtS  is  the  specific  heat 
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for  species  s.  The  total  differentials  are  evaluated  as  func¬ 
tion  of  the  conserved  variables  Qk,  thus  yielding: 


dp  =  (7  fe-1) 
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where  we  have  defined  Rs  =  fcs/ms  and 
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7fc-l=-^- -  (6) 

E  Pscv,s 

s=£e 

We  can  also  extract  the  frozen  speed  of  sound, 
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where  Se  =  pse.  Since  the  excited  levels  are  convected 
as  individual  species,  there  is  no  electronic  contribution 
to  the  internal  energy,  and  only  the  translational  degrees 
of  freedom  need  be  convected  with  the  total  energy,  viz. 

eh  =  *ThJ2Rs  (8) 

s^e 


such  that  jh  —  7e  =  5/3. 

In  the  current  work,  transient  effects  are  neglected  and 
Eq.  (4)  can  be  simplified  to 


^=A~'n 

ax 


(9) 


An  implicit  discretization  can  be  obtained  via  a  Taylor 
series  expansion  in  x  of  the  RHS  such  that 


34-in \ 

I -Ax  gQ  \AQ  =  AxA~1n  (10) 


in  which  the  Jacobian  is  expanded  via 

(dA-1!/) 


J  = 


dQ 
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(ii) 


For  greater  accuracy  and  stability,  the  Jacobian  has  been 
evaluated  analytically.  The  inverse  Jacobian  is  easily 
found  from  its  diagonalization  A-1  =  RA~1L,  along  with 
dA/dQ. 


III.  COLLISION  AL-R  ADI  ATI  VE  MODEL 

In  this  section,  we  describe  the  collisional-radiative 
model,  for  which  rates  are  computed  from  collision  cross 
sections  of  the  associated  processes  under  the  assump¬ 
tion  of  a  Maxwellian  electron  energy  distribution  func¬ 
tion  (EEDF).  The  CR  model  used  here  draws  upon  the 
previous  work  of  Vlcek19  and  Bultel  et  al20  both  of  which 
are  specific  to  argon.  The  model  includes  the  first  31  ex¬ 
cited  levels  of  neutral  argon,  as  well  as  the  ground  states 
of  the  singly-charged  ion.  The  population  numbers  of 
the  states  beyond  those  considered  here  are  negligible 
for  the  plasma  conditions  considered,  as  verified  by  equi¬ 
librium  considerations.  There  is  no  grouping  of  levels, 
as  the  approach  here  is  to  look  in  detail  at  the  influ¬ 
ence  of  the  electronic  states,  and  we  presently  neglect 
the  influence  of  molecular  Ar£.  As  the  CR  model  is 
being  validated  on  the  UTIAS  shock  tube  experiments, 
any  assumptions  made  herein  are  based  on  the  plasma 
conditions  of  those  experiments:  this  includes  number 
densities  up  to  1024  [m-3]  and  temperatures  just  above 
1  eV. 


A.  Processes 

The  processes  to  be  included  can  be  represented  by  the 
following  equations, 


Ar(i)  +  e  Ar(j)  +  e  (12) 

Fji 

Ar(i)  +  Ar(l)  Ar(j)  +  Ar(l)  (13) 

Lji 

Ar(i)  +  e~  Ar+  +  e~  +  e~  (14) 

Oi 

Ar(i )  +  Ar(l)  Ar+  +  e~  +  Ar(l)  (15) 

Wi 

Ar(i)  +  hvij  Ar(j)  (16) 

A  jiAji 

Arii)  +  hv  Ar+  +  e~  (17) 

A  iRi 


where  the  rate  coefficients  are  defined  by  Vlcek  19 .  We 
adopt  the  convention  j  >  i  such  that  C,j  and  Kij  repre¬ 
sent  the  rates  for  excitation  from  level  i  to  j  while  Fji , 
Lji,  and  A ji  are  the  rates  for  de-excitation  from  level 
j  to  i.  The  levels  considered  in  the  present  model  in¬ 
clude  only  those  up  to  the  3p53d  and  3p55s  manifolds 
(cf.  Table  II).  This  implies  that  ionization  and  recombi¬ 
nation  should  proceed  from  and  to  only  these  31  levels; 
although  levels  beyond  this  manifold  are  more  that  1  eV 
away  from  the  ionization  limit,  the  combination  of  small 
energy  gap  and  large  cross-section  makes  the  ionization 
from  these  higher  levels  extremely  rapid,  certainly  occur¬ 
ring  with  time  scales  much  below  the  resolution  needed 
in  the  current  research.  Thus  the  31  lower-lying  states 
considered  here  contain  the  bottleneck  to  the  ionization 
regime  and  most  of  the  energy  into  bound  states.  Note 
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that  neutral  argon  has  two  ionization  potentials  owing 
to  the  fact  that  the  levels  of  argon  are  split  according  to 
two  possible  configurations  for  the  core  angular  momen¬ 
tum  jc.  As  ionization  and  recombination  can  proceed 
only  between  levels  having  the  same  value  of  jc ,  the  two 
lowest  lying  levels  of  Ar+  must  be  accounted  for  individ¬ 
ually.  Thus  in  the  current  model,  these  two  levels  are 
treated  as  separate  species  and  convected  as  such.  Tran¬ 
sitions  between  states  with  different  jc  are  considered  as 
well,  despite  being  suppressed  relative  to  identical  core 
transitions. 

The  forward  rate  coefficients  are  computed  assum¬ 
ing  a  Maxwellian  electron  energy  distribution  function 
(EEDF)  according  to 

v 

k  =  (fcBr)2  J  a(£)£e~e/klsT  de  (18) 

where  e0  is  the  threshold  energy  and  v  is  the  mean  ther¬ 
mal  velocity, 


v  = 


(19) 


with  /i  the  reduced  mass.  The  backward  rates  are  then 
computed  from  the  principle  of  detailed  balance.  For 
excitation  processes  this  takes  the  form 


h. dex  _  i„exc  9i  ea/kBT 
ji  A/ in  c 

9j 


(20) 


where  =  Ej  —  is  the  energy  difference  between  the 
upper  and  lower  levels.  The  temperature  in  Eq.  (20)  is 
that  of  the  impacting  species.  For  ionization  and  recom¬ 
bination  processes,  the  reverse  rates  are  also  computed 
via  detailed  balance,  leading  to: 


l,rec  _ 

rl/A  -  l\lA 


-h 


(21) 


where  Ze  is  the  electron  partition  function 


Ze  =  2 


/  2irmekBTe 

V  h2 


3/2 


(22) 


Since  the  electronic  states  are  considered  separately,  the 
ratio  of  partition  functions  for  Ar  and  Ar+  is  simply 


ji  ^  9i  Ij/kyT 
Z+  ~  g+ 


(23) 


with  li  being  the  ionization  potential  of  the  ith  excited 
state,  allowing  Eq.  (21)  to  be  written  as 


urec  _  1,1 

ry*4  —  rvj 


t  J/y  1 

9+  2 


ti 2 


2'KmPkBTP 


3/2 


Ji/kBT 


(24) 


Note  that  the  exponential  temperature  dependence  in 
Eq.  (24)  is  a  function  of  the  third  body — T  is  taken  to 
be  either  the  electron  or  heavy  particle  temperature  for 


cases  in  which  the  third  body  is  an  electron  or  atom,  re¬ 
spectively.  We  can  also  use  the  Boltzmann  equilibrium 
conditions  to  define  an  ’’excitation”  temperature  which 
describes  the  actual  ratio  of  the  population  densities,  i.e.: 


T  . .  =  _ 

X,IJ  — 


(  ni/9j\ 

\ni/9i) 


-1 


(25) 


When  the  lower  level  is  the  ground  state,  we  can  simply 
denote  this  as  Txj .  Similarly,  we  can  define  a  ’’Saha  tem¬ 
perature”  to  describe  the  actual  ratio  of  ionized  species 
versus  the  neutral  population,  which  involves  solving  for 
Tg.i  the  equation: 


n+ne  9j_  =  / 2nmekBTsti\  3/2  e_/i/fcBTs,i _ 

rii  g+  V  h2  ) 

In  addition  to  species  production  rates,  the  second  mo¬ 
ment  with  respect  to  energy  is  also  necessary  to  deter¬ 
mine  the  energy  production  rate  for  photorecombination 
(radiative  capture)  processes, 

-  pOO 

R>  =  (fcEr)2  J  a(£)£2e~s/kBT  de ■  (27) 

In  practice,  all  rate  coefficients  and  their  derivatives  are 
computed  a  priori  and  tabulated  as  a  function  of  tem¬ 
perature. 


B.  Atom  impact  processes 


Neglecting  precursor  photo-ionization,  atom-atom  col¬ 
lisions  are  the  only  means  for  initial  electron  production 
behind  the  viscous  shock  and  yield  the  priming  electrons 
that  eventually  trigger  inelastic  electron  collisions.  Exci¬ 
tation  from  ground  state  in  particular  plays  a  significant 
role  in  determining  the  overall  induction  length,  i.e.  the 
distance  between  shock-front  and  electron  avalanche.  In 
addition,  the  excitation  rates  have  a  profound  influence 
on  the  dynamics  of  the  simulation.  It  will  be  shown  in 
the  second  article  that  pressure  waves  initiated  at  the 
electron  avalanche  are  exacerbated  for  undervalued  cross 
sections  while  excessive  values  can  lead  to  over  damping 
of  the  system.  The  sharp  sensitivity  of  the  convection- 
kinetics  coupling  highlights  the  need  for  accurate  deter¬ 
mination  of  the  heavy  particle  cross  sections.  We  start 
with  Drawin’s  formula21 


<(e)  =  ^a20 


LH 


kO-Ar  j-2  £ 
“S  Jij 


mH 


2  me 


e/Ii  -  1 


VflAr 


i  + 


2m  e 


mAr+m, 


i£/h~  i))' 


(28) 


where  Ih  is  the  ionization  potential  of  the  hydrogen  atom 
and  £  =  6  is  the  number  of  optical  electrons  of  argon.  For 
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i 

e(j)  [eV] 

9i 

jc 

n£[K]j 

i 

e{i)  [eV] 

9i 

Jc 

ni[K]j 

i 

0 

1 

1.5 

[Mg]3pb 

18 

13.903 

5 

1.5 

3d[3/2]2 

2 

11.548 

5 

1.5 

4s[3/2]2 

19 

13.979 

9 

1.5 

3d[7/2]4 

3 

11.624 

3 

1.5 

4s  [3/2]  i 

20 

14.013 

7 

1.5 

3d[7/2]3 

4 

11.723 

1 

0.5 

4s'[l/2]0 

21 

14.063 

5 

1.5 

3d[5/2]2 

5 

11.828 

3 

0.5 

4s' [1/2]  i 

22 

14.068 

5 

1.5 

5s  [3/2]  2 

6 

12.907 

3 

1.5 

4p[l/2]i 

23 

14.090 

3 

1.5 

5s[3/2]i 

7 

13.076 

7 

1.5 

4p[5/2]3 

24 

14.099 

7 

1.5 

3d[5/2]3 

8 

13.095 

5 

1.5 

4p[5/2]2 

25 

14.153 

3 

1.5 

3d[3/2]i 
3d' [5/2]  2 
3d' [3/2]  2 
3d' [5/2]  3 
5s'[l/2]0 
5s' [1/2]  i 
3d' [3/2]  i 

9 

13.153 

3 

1.5 

4p[3/2]i 

26 

14.214 

5 

0.5 

10 

13.172 

5 

1.5 

4p[3/2]2 

27 

14.234 

5 

0.5 

11 

13.273 

1 

1.5 

4p[l/2]o 

4p'  [3/2]  ]. 
4p'[3/2]a 
4p'[l/2]i 

28 

14.236 

7 

0.5 

12 

13.283 

3 

0.5 

29 

14.241 

1 

0.5 

13 

13.302 

5 

0.5 

30 

14.255 

3 

0.5 

14 

13.328 

3 

0.5 

31 

14.304 

3 

0.5 

15 

13.480 

1 

0.5 

4p'[l/2]o 

oo 

15.760 

4 

1.5 

[Mg]3p5 

16 

13.845 

1 

1.5 

3d[l/2]0 

oo' 

15.937 

2 

0.5 

[Mg]3p5 

17 

13.864 

3 

1.5 

3d[l/2]i 

TABLE  II.  Argon  levels  considered  in  current  CR  model. 


The  linear  /3-parameters  are  provided  in  Tables  III  and 
IV.  Note  the  relative  strength  of  A jc  =  0  transitions 
compared  to  the  4s  —  4s'  (A jc  ^  0)  transition. 

Less  sensitive  are  the  ionization  cross  sections — we  re¬ 
mark  here  that  decreasing  these  values  by  an  order  of 
magnitude  did  not  have  any  significant  effect  on  the  in¬ 
duction  length  or  flow  dynamics.  The  cross  section  for 
ionization  from  ground  state  argon  is  taken  from22, 

<(e)  =  1.8  x  10”25(e  -  15.760)1'3  [m2],  (32) 


i 

j 

Pii  [m2/eV] 

fij 

i 

3 

~ Tioincr23- 

6.78  x  10~U2 

i 

5 

3.95  x  10“23 

2.56  x  10~01 

i 

17 

9.58  x  10-26 

1.00  x  10^03 

i 

23 

3.10  x  10“24 

3.40  x  10~°2 

i 

25 

8.55  x  10”24 

9.51  x  10^°2 

i 

30 

1.59  x  10“24 

1.80  x  10”°2 

i 

31 

6.92  x  10-24 

7.94  x  10~°2 

TABLE  III.  Atom  impact  excitation  parameters  for  allowed 
transitions  from  ground  state  Ar. 


i 

j 

Pti  [m2/eV] 

2 

3 

1.79  x  10'24 

2 

4 

4.80  x  10~26 

2 

5 

4.80  x  10~26 

3 

4 

4.80  x  10'26 

3 

5 

4.80  x  10-26 

4 

5 

1.79  x  10'24 

TABLE  IV.  Atom  impact  excitation  parameters  for  intra-4s 
transitions. 


while  the  formula  of  Drawin  has  been  applied  for  all  other 
levels  (for  further  details,  the  reader  is  directed  to20  and19 
and  the  references  therein). 


of(e)  =  47r°o 


/  Ih  \  ^ 2 

V  h  )  mH 


2  mP 


e/Ii  -  1 


TflAr 


me 


i  + 


2  me 


mAr+m, 


ie/Ii-  1))' 


(33) 


the  energy  ranges  under  consideration,  Eq.  (28)  is  well- 
approximated  by  a  linear  function, 


where 


of,-  0)  =  Pij(£~£ij), 


q*  _A _ 2{Iii)2/-2f  2me 

h j  —47 raQ  3  £  Jij 

J  e\-  mH 


(29) 

(30) 


The  inner  3p54s  manifold  transitions,  however,  take  the 
form20 


<(£)  =  Pii 


£  £{j 

c-2.26 


(31) 


C.  Electron  impact  processes 


Once  a  sufficient  number  of  priming  electrons  have 
been  generated  by  heavy  particle  impact  and  their 
temperature  increased  by  thermalization,  electron  im¬ 
pact  processes  begin  to  become  important  and  even¬ 
tually  dominate  the  kinetics.  The  cross  sections  for 
electron-impact  excitation  are  due  to  Zatsarinny  and 
Bartschat23,24.  These  include  excitation  from  the  ground 
and  4s  levels  to  all  levels  below  the  5 p  manifold.  For  all 
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parameter 
for  Eq.  (36) 

s 

Valence  electron  shell 
p  d 

/ 

a 

1.06 

2 

3/2 

3/2 

b 

0.23 

1 

3 

1 

c 

1 

1 

2/3 

2/3 

d 

1 

1 

1 

1 

TABLE  V.  Electron-impact  ionization  parameters  as  taken 
from26. 


Valence  electron 
shell,  ni 

rni  [A] 

9n£  X  In£  [gV] 

4  s 

2.49 

7.40 

5  s 

6.35 

6.35 

4 p 

3.40 

31.00 

3  d 

4.36 

13.60 

TABLE  VI.  Radii  of  Ar  valence  electron  and  reduced  weight¬ 
ing  factors  for  £  =  1  as  taken  from26  and25. 


other  transitions,  Drawin’s  formulas 


with  the  Debye  length  (neglecting  the  effect  of  ions), 


feQkBTe\1/2 
1  Nee2  ) 


(38) 


For  the  shock  tube  conditions  under  consideration,  the 
lowering  of  the  ionization  potential  was  found  to  be  rel¬ 
atively  small  (0.08),  yielding  a  drop  of  0.08  eV  for  an 
electron  temperature  of  1  eV  and  number  density  of 
1023  m  3.  As  a  result,  the  effects  of  pressure  ionization 
were  neglected. 


D.  Photorecombination 

In  the  absence  of  a  third  body,  the  energy  released  in 
recombination  is  liberated  as  radiation.  Photorecombi¬ 
nation  is  a  significant  loss  mechanism  and  plays  an  im¬ 
portant  role  in  radiative  cooling.  The  cross  section  for 
photorecombination  can  be  found  from  the  cross  section 
for  photoionization  under  equilibrium  conditions28, 


2  ( uij  1) 


fa  (ft)" ) 


crf-(£)  =  47ra; 


(34) 

have  been  used  systematically  with  =  1,  erf)  =  1, 


and  erf)  =  1. 


Cross  sections  for  electron-impact  ionization  from  the 
excited  levels  have  been  determined  based  on  the  work 
of  Deutsch  et  al.25: 


o-i(e)  =  gniTrr2n£ntf{,e)  (35) 


<7c,i(^) 


9n  h2v2 
g+  m2v2c 2 


(39) 


Utilizing  the  relation  hv  =  mev2 /2+Ei  =  e+£j,  the  cross 
section  associated  with  the  ground  state  is  given  by 


eri(e) 


9i  (e  +  £i)2 
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2.8  x  IQ-20 
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(40) 


where  9nt  are  the  reduced  weighting  factors,  rn£  are  the  while  all  others  are  computed  using 
radii  of  the  valence  electron,  and 


fi(e)  =dJ 


e/h  -  1 


b  +  c 


e/ h  +  1 

(l-^in(2.7  +  (e/Ji-l)V2) 


(36) 


The  necessary  parameters  for  Eqs.  (35)  and  (36)  are  given 
in  Tables  V  and  VI. 

Particular  attention  should  be  paid  when  considering 
ionization  and  recombination  in  plasmas  defined  by  rel¬ 
atively  high  densities,  as  non-ideal  effects  for  which  in¬ 
teractions  between  particles  can  lead  to  an  effective  low¬ 
ering  of  the  ionization  potential.  As  a  consequence,  the 
rates  must  be  considered  as  functions  of  pressure  as  well. 
Griem2'  used  Debye  theory  to  predict  a  decrease  in  the 
ionization  potential  that  is  inversely  proportional  to  the 
Debye  length, 


(to  +  l)e2 
47re0A_D 


(37) 
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gi  (e  +  Ej)2 

3+  2  emec2 


1CT22 
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5/2 
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E  +  Ei 


0  <  e  <  0.59 Ih  —  Ei 
e  >  0.59 Ih  —  Ei 


(41) 


The  parameter  7 i(npqn^)  takes  the  values  0.0763,  0.0458, 
0.0305,  and  0.0915  for  i  =  2,  3,  4,  and  5,  respectively19. 

Based  on  nominal  shock  tube  plasma  conditions,  an 
estimation  of  the  mean-free  path  for  photoionization  from 
the  ground  state  is  on  the  order  of  100  /im.  In  the  case 
of  the  4s  levels,  their  lower  populations  combined  with 
reduced  cross  sections  result  in  mean-free  paths  on  the 
order  of  1  m.  Consequently,  free-bound  radiation  into  the 
ground  state  is  assumed  to  be  locally  absorbed,  while  the 
plasma  is  assumed  to  be  optically  thin  to  transitions  from 
the  excited  levels. 
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lower  upper  mean-free 

level  level  path  [m] 


4s  [3/2]  1 

3.59  xlO"2 

■ip6 

4s' [1/2]! 

9.47  xlO-3 

■ip6 

3d[l/2]i 

2.41 

■ip6 

5s[3/2]i 

7.11  xlO-2 

■ip6 

3d[3/2]i 

2.54  xlO-2 

■ip6 

5s' [1/2]  1 

1.34  xlO-1 

■ip6 

3d'  [3/2]  1 

3.05  xlO"2 

TABLE  VII.  Estimated  mean-free  paths  for  bound-bound 
transitions  to  ground  state  for  T  =  1  eV  and  n3pe  =  1024  m"3. 


F.  Bremsstrahlung  emission 


Free-free  transitions  have  been  incorporated  in  the 
model  via  Kramer’s  formula28  for  Bremsstrahlung  emis¬ 
sion.  Thus  the  corresponding  term  in  the  vector  Q  of 
Eq.  1  is: 


dEe  1  16-7T2  veZ^se6g 

dt  J0R  3^  me/i(47re0c)3 

=  -1.42  x  10"40Z2ffTe1/2n+ne 


(47) 

[J/m3  •  s] 


E.  Bound-bound  transitions 


Bound-bound  transitions  are  another  significant  source 
of  radiative  cooling  of  the  plasma.  The  photo-absorption 
cross  section  for  a  given  bound-bound  transition  may  be 
estimated  by28, 


aij  ~~ 


-  [  dv  =  2.65  x  10"6^ 


where  the  absorption  oscillator  strength  is  given  by 


T*abs 

Jij 


9j  Ajj 

9i  37  ’ 


(43) 


Contributions  to  the  parameter  7  have  been  assumed  to 
result  from  a  combination  of  natural  and  pressure  line 
broadening, 


7  7nat  4“  7col 


(44) 


given  by 


87r2e2^2 


7nat  =  q  3  =  2.47  x  10"2V  [s"1]  (45) 


and 


7coi  =  - =  2  avn 

7col 


(46) 


respectively.  The  collision  time  rcoi  and  the  collision  cross 
section  cr  have  been  tallied  for  all  argon-argon  collisions. 
The  radiation  cross  sections  Eq.  (42)  have  been  estimated 
using  the  oscillator  strengths  and  transition  probabilities 
from29  along  with  representative  plasma  conditions  of  the 
shock  tube  experiments.  The  resulting  mean-free  paths 
for  bound-bound  transitions  to  ground  state  given  in  Ta¬ 
ble  VII  are  below  the  length  scales  of  the  shock  tube 
dimensions  and  induction  length  for  most  of  the  levels. 
It  has  therefore  been  assumed  that  all  bound-bound  ra¬ 
diation  to  the  ground  level  from  the  excited  states  is  ab¬ 
sorbed.  For  all  other  transitions,  the  mean-free  path  is 
several  orders  of  magnitude  greater  than  the  dimensions 
of  the  shock  tube  and  the  associated  radiation  is  assumed 
to  escape.  This  simplified  treatment  of  the  radiation  will 
be  re-examined  in  a  future  study  and  compared  with  a 
radiation  transport  model. 


where  g  is  the  gaunt  factor  (here  taken  to  be  unity) 
and  the  effective  charge  Z2ff  is  taken  to  be  1.6711. 
Bremsstrahlung  emission  resulting  from  neutral  atoms 
is  1-2  orders  of  magnitude  less  than  for  ions  and  has 
therefore  been  neglected28.  The  plasma  is  assumed  to  be 
optically  thin  to  all  Bremsstrahlung  emission. 


G.  Elastic  collisions 


Elastic  collisions  are  incorporated  into  the  CR  implicit 
solver  as  well.  This  strong  coupling  permits  more  accu¬ 
rate  and  stable  calculations.  The  energy  transfer  between 
electrons  and  heavy  particles  is  computed  from 


=  nenn^kB(Th  -  Te)6keh  (48) 

CR  ^ 

where  6  =  2 me/mAr  and  keh  =  kei,  ken  for  Coulomb  and 
elastic  collisions  with  neutral  argon  respectively.  For  the 
latter,  the  theoretical  cross  section  of  McEachran  and 
Stauffer30  have  been  utilized  which  are  reproduced  in 
Table  VIII.  The  Coulomb  collision  rates,  kei,  have  been 
computed  using  the  energy-averaged  properties31, 
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kf'i  VpGp 


dei  =  5.85  x  10"10  —  [m2l 


Jn  A 

~rf 


(49) 

(50) 


where  aei  is  the  momentum-averaged  cross  section,  and 


/  T:i  \  1//2 

A  =  1.24  x  107  (51) 


is  the  Coulomb  logarithm. 


e  [eV] 

ae„  [1020  m2] 

£  [eV] 

ae„  [1020  m2] 

0.01 

4.0231 

0.41 

0.3396 

0.03 

2.3425 

0.51 

0.5349 

0.05 

1.5247 

0.61 

0.7381 

0.07 

1.0261 

0.71 

0.9318 

0.09 

0.6998 

0.81 

1.1139 

0.13 

0.3291 

0.91 

1.2835 

0.17 

0.1599 

1.00 

1.4280 

0.19 

0.1198 

1.50 

2.1438 

0.21 

0.0983 

2.00 

2.8318 

0.23 

0.0917 

3.00 

4.4046 

0.25 

0.0969 

4.00 

6.3502 

0.29 

0.1316 

5.00 

8.6317 

0.31 

0.1573 

7.50 

14.6559 

0.32 

0.1724 

10.00 

17.8325 

TABLE  VIII.  Momentum  transfer  cross  sections  for  electron- 
neutral  collisions30. 


H.  Rate  equations 


as  well  as  that  for  the  electron  energy, 
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These  equations  are  solved  simultaneously  with  a 
backward-Euler  method,  i.e.  are  implicit  in  time. 


The  rate  equations  for  the  ground  and  excited  states 
of  neutral  argon  are  given  by 
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(52) 


while  those  for  the  two  ground  states  of  singly-ionized 
argon  take  the  form 
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DT 
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(53) 


The  energy  rate  equations  are  most  easily  split  into  one 
for  the  heavy-particle  energy, 
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IV.  VALIDATION/CALIBRATION 
A.  Atom  Impact  Calibration 

Using  dual-wavelength  interferometry,  the  UTIAS 
experiments11  were  able  to  provide  detailed  insight  of 
the  shock  structure  in  terms  of  the  total  mass  density 
of  the  plasma  as  well  as  the  electron  number  density. 
From  this  it  was  possible  to  obtain  not  only  the  relaxation 
length  £  and  the  peak  ionization  fraction  a,  but  also  the 
slopes  of  the  plasma  properties  in  the  electron  avalanche 
and  in  the  radiative-cooling  region,  which  could  be  used 
to  validate  the  CR  model.  Experimental  conditions  for 
the  seven  UTIAS  cases  studied  here  are  provided  in  Ta¬ 
ble  IX.  Digitization  of  the  interferometry  snapshots  have 
yielded  results  of  the  relaxing  shock  which  are  presented 
in  terms  of  total  mass  density  p,  electron  number  density 
ne,  and  ionization  fraction  a  in  the  provided  references. 
As  there  exists  a  fair  amount  of  uncertainty  associated 
with  the  atom-atom  excitation  cross  sections  cra,  these 
were  varied  until  an  acceptable  agreement  between  the 
theoretically-predicted  and  experimentally-observed  re¬ 
laxation  lengths  was  achieved.  Exact  agreement  for  all 
cases  could  not  be  achieved  with  a  fixed  value  for  era. 
However,  good  agreement  was  obtained  for  the  interme¬ 
diate  cases,  corresponding  to  Mach  numbers  14.7,  15.9, 
and  16.1,  the  results  of  which  are  summarized  in  Table 
X. 

As  mentioned  in  Section  IIIB,  the  atom-atom  impact 
excitation  cross  sections  have  been  approximated  as  lin¬ 
ear  functions  of  energy  as  based  on  the  model  of  Drawin. 
Tuning  these  cross  sections  was  therefore  a  matter  of 
modifying  the  slopes  of  the  parameter  /?.  As  a  result, 
the  values  of  3*j  in  Eq.  29  were  reduced  by  a  factor  of 
25,  and  are  summarized  in  Table  XI.  In  comparison  with 
values  previously  obtained  in  the  literature  (cf.  Table  I), 
it  is  evident  that  the  cross  sections  used  here  are  at  the 
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Case 

Mas 

Po  [torr] 

To  [K] 

Ref. 

1 

15.9 

5.14 

293.6 

Tang  (1977) 

2 

16.1 

5.15 

295.9 

Tang  (1977) 

3 

16.5 

5.12 

296.6 

Brimelow  (1974) 

4 

13.0 

5.01 

296.6 

Whitten  (1977) 

5 

17.2 

2.06 

297.8 

Bristow  (1971) 

6 

13.6 

5.12 

297.0 

Brimelow  (1974) 

7 

14.7 

4.08 

297.8 

Glass  &  Liu  (1978) 

TABLE  IX.  Shock  tube  cases  studied  at  UTIAS. 


Experiment  Theory 


Case 

it  [cm] 

a* 

it  [cm] 

a* 

1 

2.00 

0.140 

2.00 

0.133 

2 

1.90 

0.150 

1.80 

0.141 

3 

1.80 

0.160 

1.53 

0.154 

4 

8.90 

0.060 

10.9 

0.058 

5 

1.60 

0.180 

2.31 

0.180 

6 

8.00 

0.074 

6.96 

0.072 

7 

4.40 

0.106 

4.35 

0.103 

TABLE  X.  Comparison  between  experimental  and  steady- 
state  numerical  relaxation  lengths  for  ionizing  shock  in  argon. 


<M 
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lower  end  of  the  spectrum,  matching  most  closely  with 
the  values  obtained  by  McLaren  and  Hobson6. 


FIG.  1.  Ionizing  shock  structure  as  detailed  by  (a)  p  and  ne, 
and  (b)  a  for  case  1:  pa  =  5.14  torr,  T0  =  293.6  K,  Ma  =  15.9. 


B.  Relaxation  Region 

The  results  using  these  cross  sections  for  case  1  given 
in  Fig.  1  show  excellent  agreement  with  the  experimen¬ 
tal  results.  In  an  effort  to  quantify  the  influence  of  the 
upper  levels  on  the  overall  shock  structure,  simulations 
were  first  performed  taking  into  account  the  4s  mani¬ 
fold  of  Ar  exclusively,  followed  by  systematic  inclusion  of 
higher  manifolds,  including  the  4 p,  5s  and  3 d.  Numerical 
results  obtained  with  levels  of  only  the  4s  manifold  given 
in  Figure  1  show  satisfactory  prediction  of  the  induction 
length,  but  poor  reproduction  of  the  radiative  cooling 
region  as  indicated  by  the  slow  drop-off  in  electron  num¬ 
ber  density,  indicating  an  under-prediction  in  radiative 
losses.  This  is  as  expected  as  the  plasma  has  been  as- 


i 

3 

P*j  [m2/eV] 

i 

3 

9.35  x  10~25 

i 

5 

3.36  x  10~24 

i 

17 

8.14  x  10~27 

i 

23 

2.64  x  10~25 

i 

25 

7.27  x  10~25 

i 

30 

1.35  x  10~25 

i 

31 

5.88  x  10~25 

TABLE  XI.  Tuned  atom  impact  excitation  parameters  for 
allowed  transitions  from  ground  state  Ar. 


sumed  to  be  optically-thick  for  all  transitions  to  ground 
state  Ar  (cf.  Sections  HID,  HIE),  resulting  in  no  radia¬ 
tive  losses  due  to  bound-bound  transitions  since  intra-4s 
transitions  are  forbidden.  As  a  result,  all  radiative  losses 
in  this  case  were  due  to  free-free  and  free-bound  transi¬ 
tions  only.  Also  apparent  is  the  effect  of  the  upper  levels 
on  the  induction  length  -addition  of  the  upper  levels  pro¬ 
motes  ladder  climbing  in  the  induction  zone,  facilitating 
ionization  and  onset  of  the  electron  avalanche.  It  should 
be  noted  that  the  ionization  fractions  as  derived  from  the 
experimental  data  for  this  particular  case,  Ma  =  15.9,  as 
well  as  for  the  case  of  Ma  =  16.1  in11  do  not  coincide  with 
the  respective  total  mass  and  electron  number  densities 
provided.  As  a  consequence,  the  experimental  ionization 
fractions  plotted  for  comparison  have  been  obtained  di¬ 
rectly  from  the  plasma  densities. 

Details  of  the  electron  and  heavy-particle  temperature 
profiles  along  with  the  Boltzmann  and  Saha  equivalence 
temperatures  provided  in  Figure  2  help  illustrate  the  in¬ 
fluence  of  the  various  CR  processes  and  to  separate  the 
shock  structure  into  five  distinct  regions.  The  first  region, 
easily  identifiable  by  a  sharp  spike  in  electron  tempera¬ 
ture  just  behind  the  compression  shock,  indicates  the  ini¬ 
tial  production  of  priming  electrons.  The  spike  is  quickly 
followed  by  a  dip  in  electron  temperature  (region  II) ,  sig¬ 
naling  a  shift  from  heavy-particle  dominated  kinetics  to 
electron- impact  processes.  The  excitation  temperatures 
based  on  a  Boltzmann  distribution  (see  Eq.  (25))  and 
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FIG.  2.  Plot  of  heavy-particle  and  electron  temperatures 
along  with  excitation  temperatures  of  excited  states  for  Mach 
15.9  shock. 


FIG.  3.  Net  power  transfer  due  to  (a)  elastic  and  (b)  inelastic 
collisions. 

measured  with  respect  to  the  ground  state,  surpass  the 
electron  temperature,  indicating  over-population  of  the 
excited  states  and  confirming  that  the  priming  electrons 
are  generated  via  Ar-Ar  inelastic  collisions.  The  evolu¬ 
tion  of  Te  presumably  results  from  a  competition  between 
elastic  exchanges  with  the  hot  atoms  and  energy  loss  due 
to  excitation  and  ionization.  To  confirm  this  interpreta¬ 
tion,  we  monitored  the  effective  rates  of  change  to  the 
electron  energy  and  density  due  to  both  elastic  and  in¬ 
elastic  processes,  shown  in  Figure  3.  It  is  clear  that  while 
elastic  collisions  gradually  increase  the  electron  energy, 


FIG.  4.  Accumulated  net  electron  number  density  production 
separated  by  manifold. 


inelastic  collisions  act  as  a  significant  energy  sink  term. 
The  point  at  which  the  electron  temperature  starts  to  dip 
(start  of  region  II)  corresponds  to  that  for  which  electron 
impact  ionization  from  the  upper  levels  starts  to  dom¬ 
inate  over  the  heavy-particle  impact  process  (compare 
C,F  and  V,W  curves  in  Fig.  3).  After  this  small  drop, 
the  electron  temperature  recovers  and  increases  steadily 
as  energy  is  transferred  from  the  heavy-particles  to  the 
free  electrons  at  a  slightly  higher  rate  than  which  the 
electrons  lose  energy  through  inelastic  collisions.  In  this 
region  (III),  electron  production  begins  to  ramp  up  expo¬ 
nentially,  triggering  rapid  thermalization  with  the  heavy- 
particles.  Note  that  the  excitation  temperature  declines 
as  well,  and  stabilizes  significantly  below  the  electron 
temperature.  This  can  be  explained  by  rapid  ionization 
from  the  excited  levels,  causing  under-population  of  the 
excited  levels,  and  hence  the  deviation  from  Te.  It  is  also 
worth  pointing  out  that  the  start  of  this  region  approxi¬ 
mately  coincides  with  the  dominance  of  electron-impact 
ionization  (compare  S,0  and  V,W  curves)  and  soon  af¬ 
ter  the  dominance  of  electron-ion  thermalization.  Thus, 
as  more  electrons  (and  ions)  are  generated  by  electron 
impact  ionization,  the  rate  of  thermalization  increases, 
which  further  enhances  the  rate  of  electron  production. 
This  highly  non-linear  process  inevitably  results  in  the 
electron  avalanche  (region  IV),  when  sufficient  electrons 
are  produced  to  affect  the  temperature  of  the  neutral 
bath  itself,  which  rapidly  decays.  The  avalanche  results 
in  complete  thermalization  between  the  heavy-particles 
and  electrons  (region  V),  as  the  energy  reserves  of  the 
heavy  particles  becomes  exhausted.  All  levels  finally 
reach  Boltzmann  and  Saha  equilibrium  and  at  that  point 
the  plasma  is  well-described  by  a  single  temperature. 

Note  that  the  spread  of  excitation  temperatures 
throughout  the  shock  structure  is  quite  small,  indicating 
a  Boltzmann  equilibrium  among  the  excited  levels  with 
respect  to  the  ground  state,  due  to  a  rapid  exchange  be¬ 
tween  the  excited  levels.  Thus,  the  plasma  may  be  well 
approximated  by  a  three-temperature  model  under  such 
conditions10. 
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FIG.  5.  Individual  contributions  to  radiative  losses 


FIG.  6.  Effect  of  discretization  on  Saha  temperature:  (a) 
Ax  =  25/rm  (b)  Ax  =  2.5/xm  (c)  Ax  =  0.25/rm  . 


Since  the  evolution  of  the  electron  temperature  follows 
from  the  balance  between  thermalization  and  inelastic 
losses,  the  elastic  collisions  can  have  a  significant  effect 
on  the  overall  relaxation  length — their  inclusion  requires 
reduction  of  the  atom-impact  cross  sections  by  one  half 
to  obtain  lengths  comparable  to  cases  when  such  colli¬ 
sions  are  not  considered — and  proper  treatment  of  the 
Coulomb  collisions  is  required.  The  electron  production 
rates  shown  in  Fig.  4  corroborate  the  mechanism  of  pro¬ 
duction  of  priming  electrons  from  atom-impact  ionization 
which  is  then  quickly  eclipsed  by  electron-impact  pro¬ 
cesses.  Furthermore  we  verify  that  ionization  proceeds 
most  prevalently  from  the  upper  levels  via  the  ladder¬ 
climbing  mechanism.  Inelastic  collisional  exchange  be¬ 
tween  the  excited  levels  is  very  fast,  leading  to  a  Boltz¬ 
mann  distribution  with  a  nearly  common  Tx. 

The  radiative  power  density  is  shown  in  Fig.  5(a). 
Comparing  with  the  collisional  rates  of  Fig.  3,  it  can  be 
seen  that  the  power  is  much  smaller  in  the  induction 
zone  (<  2cm)  as  expected,  since  the  electron  density  is 
low.  At  and  behind  the  avalanche,  the  overall  radiation 
is  increased  by  several  orders  of  magnitude,  with  the  line 
transitions  dominating  the  radiative  losses.  However,  the 
strongest  lines  are  for  transitions  to  the  ground  state, 
which  are  assumed  self-absorbed  (escape  factor  A ^  =  0). 
If  we  account  for  the  self-absorption  of  both  the  bound- 
bound  and  free-bound  transitions  to  the  ground  state, 
the  net  radiative  power  would  become  much  smaller. 

Figure  6  illustrates  the  effects  on  the  avalanche  re¬ 
gion  if  the  collisional  time  scales  are  not  fully  resolved. 
For  a  coarse  discretization,  the  Saha  temperatures  of  the 
higher  states  (green,  red  and  blue  solid  curves)  overshoot 
Te  and  relaxes  with  T^,  indicating  over- ionization  in  the 
avalanche  region  w.r.t.  electron-impact  collisions.  Since 
the  electron-impact  cross  sections  are  larger  than  their 
atom-impact  counterparts,  the  Saha  temperature  should 
instead  relax  to  Te  as  in  6(c).  This  unphysical  behavior 
disappears  once  the  fast-time  scales  are  resolved  (which  is 
equivalent  to  spacing  resolution  for  the  steady-state  cal¬ 
culation).  Because  the  higher  states  have  a  low  absolute 
population  and  exchange  little  energy  with  the  contin¬ 


uum  (small  ionization  potential),  the  solution  remains 
stable  and  the  basic  shock  structure  unchanged  even  if 
under-resolved . 


C.  Electron  Heat  Conduction 


Due  to  the  high  thermal  velocities  of  the  free  electron 
gas  combined  with  large  temperature  gradients,  the  rates 
of  electron  heat  conduction  in  the  relaxation  region  can 
be  significant.  The  change  in  the  conservative  variables 
may  be  expressed  in  operator-split  form  as 


dEe  _d_(  dTA 
dt  cond  dx  V  6  dx  J  ’ 


(64) 


with  the  form  of  the  electron  thermal  conductivity  taken 


as 


32 


Ke 


5nekfrTe 
2  meve 


(65) 


where  the  total  electron  collision  frequency  includes  all 
species,  ve  =  vee  +  vei  +  i/en.  We  have  used  both  the 
energy-averaged  electron-electron  and  electron-ion  colli¬ 
sion  frequencies, 


4-\/27r  /  me  \  3^2  /  e2 
3  \kBTeJ  \4Tre0me 


In  A, 


(66) 


and  z/ee  =  2uei,  along  with  the  electron-neutral  collision 
frequency  Den  =  nhken, 

v 

ken  =  (, kBT  )2  i  a(£)£e~£/kBTed£  (67) 

Eq.  (64)  was  solved  implicitly  in  x  using  operator¬ 
splitting,  assuming  the  shock  front  to  be  adiabatic.  This 
is  consistent  with  our  assumption  that  there  are  no  pre¬ 
cursor  electrons.  At  the  opposite  end  of  the  domain, 
where  the  slope  of  the  radiative-cooling  regime  dT /dx 
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FIG.  7.  Effect  of  electron  heat  conduction  on  temperature  in 
the  relaxation  region. 


is  dominated  by  radiative  processes,  we  have  fixed  the 
slope  of  Te,  thereby  assuming  a  constant  heat  flux.  The 
effect  of  heat  conduction  applied  to  case  1  can  be  seen 
in  Fig.  7.  Most  notably,  the  sharp  gradient  in  Te  just 
behind  the  shock  has  been  eliminated  followed  by  a  large 
reduction  in  the  trough  created  by  electron-impact  pro¬ 
cesses  following  the  overshoot.  Since  the  electron  density 
is  quite  low  in  this  priming  region,  the  amount  of  en¬ 
ergy  diffused  has  a  trivial  effect  on  the  overall  relaxation 
length  and  temperature  profiles,  which  remain  essentially 
unchanged.  It  is  clear  from  comparison  of  Figures  7(b) 
and  3(a)  that  the  energy  transferred  through  conduction 
is  roughly  two  orders  of  magnitude  lower  than  the  en¬ 
ergy  trasferred  through  Coulomb  collisions.  Neverthe¬ 
less,  the  smoothing  of  the  electron  temperature  profile  in 
the  priming  region  can  be  important  from  a  numerical 
standpoint:  because  of  the  highly  nonlinear  dependence 
of  the  kinetics  on  temperature,  the  sharp  gradient  im¬ 
poses  a  strict  time-step  limitation  on  the  kinetics,  which 
translates  into  a  small  spatial  stp  (approximately  33/rm, 
equivalent  to  about  600  cells  over  the  entire  relaxation 
length  for  a  constant-spacing  grid).  This  becomes  in¬ 
creasingly  important  as  more  of  the  upper  levels  are  in¬ 
cluded  as  they  tend  to  increase  the  slope  of  Te  just  behind 
the  shock. 


D.  Other  Experiments 

Results  for  cases  2-5  are  provided  in  Fig.  8-11.  Both 
the  Mach  15.9  and  16.1  cases  are  in  excellent  agree¬ 
ment,  which  indicates  that  the  numerical  solution  of  the 
collisional-radiative  is  consistent.  However,  there  is  a 
slight  under-prediction  of  the  induction  length  for  the 
Mach  16.5  case,  and  an  over-prediction  for  the  Mach  13 
and  17.2  case.  This  non-monotonic  behavior  could  not 
simply  be  explained  from  the  collisional-radiative  model. 
There  are,  however,  several  factors  which  can  influence 
the  experimental  results.  As  noted  by  Enomoto33,  the 
boundary  layer  interaction  with  the  shock  structure  ini¬ 
tiates  premature  onset  of  the  electron  avalanche,  thus 
shortening  the  overall  relaxation  length.  Glass  and  Liu 
showed  for  a  Mach  13.6  shock,  the  boundary  layer  signif¬ 
icantly  shortened  the  relaxation  length  more  than  4  mm 
away  from  the  wall.  For  the  5  cm  diameter  ANU  shock 
tube,  this  represents  a  significant  portion  of  the  flow  that 
would  be  affected  by  such  a  boundary  layer.  Another  fac¬ 
tor  is  the  presence  of  impurities;  as  noted  also  by  Glass 
and  Liu11  and  verified  by  one-dimensional  computations 
with  a  simplified  model10,  molecular  impurities  (includ¬ 
ing  traces  of  water  vapor)  can  have  a  strong  effect  in 
reducing  the  induction  length.  Both  of  these  can  be  esti¬ 
mated  with  further  computations,  extending  the  scheme 
to  include  the  Navier-Stokes  equations  for  boundary  layer 
development,  and  molecular  chemistry.  Yet  another  fac¬ 
tor  to  consider  is  the  unsteadiness  of  the  flow;  the  exper¬ 
imental  traces  in  figures  1  and  8-11  clearly  show  signifi¬ 
cant  fluctuations.  Accounting  for  the  occasional  shift  in 
induction  length,  the  computational  profiles  of  the  mass 
density  and  electron  density  or  ionization  fraction  show  a 
remarkable  agreement  with  the  mean  of  the  experimen¬ 
tal  data.  We  will  show  in  the  companion  paper  that  the 
flow  exhibits  very  particular  unsteady  fluctuations,  which 
lead  to  the  observed  profiles,  as  well  as  effects  on  the  in¬ 
duction  length  which  can  bring  our  results  closer  to  the 
experimental  data. 

Research  on  strongly  ionizing  shocks  in  argon  were  also 
carried  out  at  the  Australian  National  University  (ANU) 
by  Houwing  et  al34,  who  was  able  to  measure  the  popu¬ 
lation  of  the  metastable  states  behind  the  shock  through 
spatially-resolved  hook  interferometry.  We  also  simu¬ 
lated  these  experiments,  choosing  similar  flow  conditions 
and  using  the  same  data-set  of  cross-sections. 

The  metastable  populations  behind  a  Mach  13.6  shock 
are  presented  in  Figure  12(a),  at  the  top.  The  cur¬ 
rent  analysis  accurately  reproduces  the  peaks  of  both 
the  4s  [3/2]  2  and  4s'  [1/2]  2  curves  as  well  as  their  slopes 
(the  rate  of  decay)  in  the  radiative  cooling  regime.  Such 
agreement  is  only  possible  when  the  4 p,  5s,  and  3d  are  in¬ 
cluded  in  the  model  in  addition  to  the  4s  levels,  otherwise 
the  population  of  the  metastable  states  is  greatly  over¬ 
predicted  as  in  Fig.  12(b), presenting  yet  another  strong 
evidence  for  the  ladder  climbing  process.  The  relaxation 
length  in  Figure  12(a),  however,  is  again  slightly  over¬ 
predicted,  even  when  experimental  uncertainties  are  ac- 
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FIG.  8.  Ionizing  shock  structure  as  detailed  by  (a)  p  and  ne, 
and  (b)  a  for  case  2:  pa  =  5.15  torr,  Ta  =  295.9  K,  Ma  =  16.1. 


x  [cm] 


FIG.  10.  Ionizing  shock  structure  as  detailed  by  (a)  p  and 
ne,  and  (b)  a  for  case  4:  pa  =  5.01  torr,  Ta  =  296.6  K,  Ma 
=  13.0. 


FIG.  9.  Ionizing  shock  structure  as  detailed  by  (a)  p  and  ne, 
and  (b)  a  for  case  3:  pa  =  5.12  torr,  T0  =  296.6  K,  Ma  =  16.5. 


x  [cm] 

FIG.  11.  Ionizing  shock  structure  as  detailed  by  (a)  p  and 
ne,  and  (b)  a  for  case  5:  p0  =  2.06  torr,  Ta  =  297.8  K,  Ma 
=  17.2. 
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FIG.  12.  Metastable  populations  behind  a  strongly  ionizing 
shock  traveling  at  4.35  ±  0.5  km/s  in  pure  argon  at  10.0  ± 
0.5  torr,  294  K. 

counted  for.  In  fact,  the  current  analysis  over-predicts 
the  relaxation  length  for  all  cases  presented  by  Houwing 
et  al.  when  using  the  atom-atom  impact  cross  sections 
calibrated  for  the  UTIAS  experiments.  The  reason  for 
such  a  consistent  error  is  likely  due  to  the  size  of  the 
shock  tube  used  in  the  experiments,  as  noted  earlier,  i.e. 
the  effect  of  developing  boundary  layers. 

V.  CONCLUSIONS 

A  detailed  collisional-radiative  model  has  been  intro¬ 
duced  which  is  able  to  accurately  reproduce  the  struc¬ 
ture  of  strongly  ionizing  shocks  in  argon.  Compari¬ 
son  between  steady-state  computations  and  experimental 
shock-tube  data  in  the  Mach  13-18  range  shows  that  the 
model  provides  a  very  good  agreement  with  the  density 
profiles  and  ionization  fraction,  and  is  able  to  reproduce 
the  population  of  the  metastable  states,  provided  that  the 
collisional-radiative  model  includes  at  least  the  5s  and  3d 
manifolds  of  neutral  argon.  We  have  also  analyzed  the  re¬ 
spective  contributions  of  the  various  elastic,  inelastic  and 
radiative  exchange  processes  and  shown  strong  evidence 
in  support  of  the  ladder  climbing  mechanism  of  ioniza¬ 
tion,  in  contrast  to  the  two-step  model  of  excitation  to 
and  ionization  from  the  4s  manifold.  Our  simulations 
indicate  that  ionization  proceeds  more  prevalently  from 
the  upper  levels,  and  the  ladder  climbing  process  is  also 


more  efficient  as  it  results  in  reduced  relaxation  lengths. 
However,  the  inelastic  collisions  between  the  upper  states 
is  sufficiently  rapid  that  a  single  excitation  temperature 
could  be  considered  a  reasonable  approximation  in  the 
induction  region.  Electron  heat  conduction  has  a  negli¬ 
gible  effect  on  the  relaxation  length,  but  can  reduce  time 
step  limitations  imposed  by  the  kinetics  in  the  electron 
priming  region.  While  we  have  been  able  to  validate  the 
model  against  shock-tube  data,  there  are  remaining  un¬ 
certainties;  these  will  be  addressed  with  further  exten¬ 
sion  of  the  model.  In  a  companion  paper,  we  address 
two-dimensional  and  unsteady  effects,  while  future  work 
will  include  boundary  layer  effects  and  more  detailed  and 
complete  kinetics,  as  well  as  radiation  transport. 
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